// WaveProb.C

#include "WaveProb.H"

void WaveProb::insert(int r, int num, double *v, int *c)
{
  if (firstFill)
    J->InsertGlobalValues(r, num, v, c);
  else
    J->ReplaceGlobalValues(r, num, v, c);
}

double WaveProb::newton(double guess)
{
  double y = guess;
  double fy = fGammaGen(y);
  if (fabs(fy) > 1e-25)
  {
    double d = -fy/dfGammaGen(y);
    while (fabs(d)/(fabs(y) + 1e-5) > 1e-8)
    {
      y += d;
      fy = fGammaGen(y);
      d = -fy/dfGammaGen(y);
    }
  }
  return y;
}

double WaveProb::fGammaGen(double lambda)
{
  double q = 0.0;
  double intGamma = 0.0;
  double hp = fabs(P0) / M;
  double next;
  for (int r=0; r<M-1; ++r)
  {
    intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
    q -= (hp/2)*pow(lambda + 2*intGamma, -1.5);
    next = intGamma - (Gamma((r+1)*hp)+ Gamma((r+2)*hp))*hp/2;
    q -= (hp/2)*pow(lambda + 2*next, -1.5);
  }
  int r = M - 1;
  intGamma -= (Gamma(r*hp) + Gamma((r+1)*hp))*hp/2;
  q -= hp*pow(lambda + 2*intGamma, -1.5);
  return g*q + 1.0;
}

double WaveProb::dfGammaGen(double lambda)
{
  double q = 0.0;
  double intGamma = 0.0;
  double hp = fabs(P0) / M;
  double next; 
  for (int r=0; r<M-1; ++r)
  {
    intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
    q -= (hp/2)*pow(lambda + 2*intGamma, -2.5);
    next = intGamma - (Gamma((r+1)*hp)+ Gamma((r+2)*hp))*hp/2;
    q -= (hp/2)*pow(lambda + 2*next, -2.5);
  }
  int r = M-1;
  intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
  q -= hp*pow(lambda + 2*intGamma, -2.5);
  return -1.5*g*q;
}

bool WaveProb::computeF(const Epetra_Vector &x, Epetra_Vector &f)
{
  double hp = fabs(P0)/M;
  double hq = 2*M_PI/N;
  int r, c;
  int knw, kn, kne, kw, k, ke, ksw, ks, kse;


  // top row
  f[0] = 1 + pow((x[1] - x[N-1])/(2*hq), 2.0) +
         (2*g*x[0] - Q)*pow((x[0]-x[N])/hp, 2.0);

  for (c=1; c<N-1; ++c)
    f[c] = 1 + pow((x[c+1] - x[c-1])/(2*hq), 2.0) +
           (2*g*x[c] - Q)*pow((x[c]-x[N+c])/hp, 2.0);

  f[N-1] = 1 + pow((x[0] - x[N-2])/(2*hq), 2.0) + 
           (2*g*x[N-1] - Q)*pow((x[N-1]-x[2*N-1])/hp, 2.0);


  // leftmost column
  for (r=1; r<M-1; ++r)
  {
    kn  = (r-1)*N;
    k   = r*N;
    ks  = (r+1)*N;
    ke  = k + 1;
    kw  = ks - 1;
    kne = kn + 1;
    knw = k - 1;
    kse = ks + 1;
    ksw = (r+2)*N  -  1;

    f[k] = (1.0 +
           pow((x[ke]-x[kw])/(2*hq), 2.0)) *
           (x[kn] - 2*x[k] + x[ks])/(hp*hp)
           -
           2*(x[kn] - x[ks])/(2*hp) *
           (x[ke] - x[kw])/(2*hq) *
           (x[kne] - x[knw] - x[kse] + x[ksw])/(4*hq*hp)
           +
           pow((x[kn] - x[ks])/(2*hp), 2.0) *
           (x[ke] - 2*x[k] + x[kw])/(hq*hq)
           +
           Gamma(hp*r) *
           pow((x[kn] - x[ks])/(2*hp), 3.0);
  }          

  
  // inner cells
  for (r=1; r<M-1; ++r)
  {
    for (c=1; c<N-1; ++c)
    {
      k   = r*N      +  c;
      ke  = k+1;
      kw  = k-1;
      kn  = (r-1)*N  +  c;
      kne = kn + 1;
      knw = kn - 1;
      ks  = (r+1)*N  +  c;
      kse = ks + 1;
      ksw = ks - 1;
      
      f[k] = (1.0 +
             pow((x[ke]-x[kw])/(2*hq), 2.0)) *
             (x[kn] - 2*x[k] + x[ks])/(hp*hp)
             -
             2*(x[kn] - x[ks])/(2*hp) *
             (x[ke] - x[kw])/(2*hq) *
             (x[kne] - x[knw] - x[kse] + x[ksw])/(4*hq*hp)
             +
             pow((x[kn] - x[ks])/(2*hp), 2.0) *
             (x[ke] - 2*x[k] + x[kw])/(hq*hq)
             +
             Gamma(hp*r) *
             pow((x[kn] - x[ks])/(2*hp), 3.0);
    }
  }


  // rightmost column
  for (r=1; r<M-1; ++r)
  {
    k   = (r+1)*N  -  1;
    kw  = k-1;
    kn  = r*N      -  1;
    knw = kn - 1;
    ks  = (r+2)*N  -  1;
    ksw = ks - 1;
    ke  = r*N;
    kne = (r-1)*N;
    kse = (r+1)*N;
    

    f[k] = (1.0 +
           pow((x[ke]-x[kw])/(2*hq), 2.0)) *
           (x[kn] - 2*x[k] + x[ks])/(hp*hp)
           -
           2*(x[kn] - x[ks])/(2*hp) *
           (x[ke] - x[kw])/(2*hq) *
           (x[kne] - x[knw] - x[kse] + x[ksw])/(4*hq*hp)
           +
           pow((x[kn] - x[ks])/(2*hp), 2.0) *
           (x[ke] - 2*x[k] + x[kw])/(hq*hq)
           +
           Gamma(hp*r) *
           pow((x[kn] - x[ks])/(2*hp), 3.0);
  }


  // bottom row
  for (c=1; c<N-1; ++c)
  {
    k   = (M-1)*N   +  c;
    kn  = (M-2)*N   +  c;
    ke  = k+1;
    kw  = k-1;
    knw = kn+1;
    kne = kn-1;
    
    f[k] = (1.0 +
           pow((x[ke]-x[kw])/(2*hq), 2.0)) *
           (x[kn] - 2*x[k])/(hp*hp)
           -
           2*(x[kn])/(2*hp) *
           (x[ke] - x[kw])/(2*hq) *
           (x[kne] - x[knw])/(4*hq*hp)
           +
           pow((x[kn])/(2*hp), 2.0) *
           (x[ke] - 2*x[k] + x[kw])/(hq*hq)
           +
           Gamma(hp*r) *
           pow(x[kn]/(2*hp), 3.0);
  }


  // bottom-left corner 
  k   = (M-1)*N;
  kn  = (M-2)*N;
  ke  = k + 1;
  kw  = k + N-1;
  kne = kn + 1;
  knw = k - 1;
  
  f[k] = (1.0 +
         pow((x[ke]-x[kw])/(2*hq), 2.0)) *
         (x[kn] - 2*x[k])/(hp*hp)
         -
         2*(x[kn])/(2*hp) *
         (x[ke] - x[kw])/(2*hq) *
         (x[kne] - x[knw])/(4*hq*hp)
         +
         pow((x[kn])/(2*hp), 2.0) *
         (x[ke] - 2*x[k] + x[kw])/(hq*hq)
         +
         Gamma(hp*r) *
         pow(x[kn]/(2*hp), 3.0);


  // bottom-right corner
  k   = M*N - 1;
  kw  = k-1;
  ke  = (M-1)*N;
  kn  = ke - 1;
  knw = kn - 1;
  kne = (M-2)*N;

  f[k] = (1.0 +
         pow((x[ke]-x[kw])/(2*hq), 2.0)) *
         (x[kn] - 2*x[k])/(hp*hp)
         -
         2*(x[kn])/(2*hp) *
         (x[ke] - x[kw])/(2*hq) *
         (x[kne] - x[knw])/(4*hq*hp)
         +
         pow((x[kn])/(2*hp), 2.0) *
         (x[ke] - 2*x[k] + x[kw])/(hq*hq)
         +
         Gamma(hp*r) *
         pow(x[kn]/(2*hp), 3.0);

  return true;
}

bool WaveProb::computeJacobian(const Epetra_Vector &x)
{
  double hp = fabs(P0)/M;
  double hq = 2*M_PI/N;
  double v;
  int k, kn, ke, ks, kw, kne, knw, kse, ksw;
  int r, c;   // index x (p-q coordinates) (e.g. r*M + c)
  int R, C;   // index J (e.g. J[R,C] )

  J->PutScalar(0.0);

  
  // D(f[0])
  R = C = 0;
  v = (2.0/(hp*hp))*(x[0]-x[N])*(3.0*g*x[0] - g*x[N] - Q); 
  insert(R, 1, &v, &C);

  C = 1;
  v = (0.5/(hq*hq))*(x[1]-x[N-1]);
  insert(R, 1, &v, &C);
  
  C = N - 1;
  v = -v;
  insert(R, 1, &v, &C);

  C = N;
  v = -(2.0/(hp*hp))*(2.0*g*x[0] - Q)*(x[0]-x[N]);
  insert(R, 1, &v, &C);


  // D(f[c]) for c=1 to N-2
  for (c=1; c<N-1; ++c)
  {
    R = C = c;
    v = (2.0/(hp*hp))*(x[c]-x[N+c])*(3.0*g*x[c]-g*x[N+c] - Q);
    insert(R, 1, &v, &C);

    C = c+1;
    v = (0.5/(hq*hq))*(x[c+1]-x[c-1]);
    insert(R, 1, &v, &C);

    C = c-1;
    v = -v;
    insert(R, 1, &v, &C);

    C = N+c;
    v = -(2.0/(hp*hp))*(2.0*g*x[c] - Q)*(x[c]-x[N+c]);
    insert(R, 1, &v, &C);
  }


  // D(f[N-1])
  R = C = N-1;
  v = (2.0/(hp*hp))*(x[N-1]-x[N+N-1])*(3.0*g*x[N-1]-g*x[N+N-1] - Q);
  insert(R, 1, &v, &C);

  C = 0;
  v = (0.5/(hq*hq))*(x[0]-x[N-2]);
  insert(R, 1, &v, &C);

  C = N-2;
  v = -v;
  insert(R, 1, &v, &C);

  C = 2*N-1;
  v = -(2.0/(hp*hp))*(2.0*g*x[N-1] - Q)*(x[N-1]-x[2*N-1]);
  insert(R, 1, &v, &C);
  


  // derivative of leftmost column
  for (r=1; r<M-1; ++r)
  {
    kn  = (r-1)*N;
    k   = r*N;
    ks  = (r+1)*N;
    ke  = k + 1;
    kw  = ks - 1;
    kne = kn + 1;
    knw = k - 1;
    kse = ks + 1;
    ksw = (r+2)*N  -  1;
    

    R = C = k;
    v = (-2/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (2/(hq*hq)*
       pow((x[kn]-x[ks])/(2*hp), 2.0));
    insert(R, 1, &v, &C);
    
    C = kw;
    v = (-0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k]+x[ks])/(hp*hp))
      +
      (1/hq)*
      ((x[kn]-x[ks])/(2*hp))*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = ke;
    v = (0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k]+x[ks])/(hp*hp))
      -
      (1/hq)*
      ((x[kn]-x[ks])/(2*hp))*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = kn;
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (0.5/(hp*hp))*
      (x[kn]-x[ks])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      +
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
    insert(R, 1, &v, &C);
    
    C = ks; 
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      +
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      -
      (0.5/(hp*hp))*
      (x[kn]-x[ks])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      -
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
    insert(R, 1, &v, &C);
    
    C = knw;
    v = (0.5/(hp*hq))*
      (x[kn]-x[ks])/(2*hp)*
      (x[ke]-x[kw])/(2*hq);
    insert(R, 1, &v, &C);
    
    C = kse;
    insert(R, 1, &v, &C);
    
    C = kne;
    v = -v;
    insert(R, 1, &v, &C);
      
    C = ksw;
    insert(R, 1, &v, &C);
    
    insert(R, 1, &v, &C);
  }


  // derivative of inner cells
  for (r=1; r<M-1; ++r)
  {
    for (c=1; c<N-1; ++c)
    {
      k   = r*N      +  c;
      ke  = k+1;
      kw  = k-1;
      kn  = (r-1)*N  +  c;
      kne = kn + 1;
      knw = kn - 1;
      ks  = (r+1)*N  +  c;
      kse = ks + 1;
      ksw = ks - 1;


      R = C = k;
      v = (-2/(hp*hp))*
          (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
          -
          (2/(hq*hq)*
          pow((x[kn]-x[ks])/(2*hp), 2.0));
      insert(R, 1, &v, &C);

      C = kw;
      v = (-0.5/(hq*hq))*
          (x[ke]-x[kw])*
          ((x[kn]-2*x[k]+x[ks])/(hp*hp))
          +
          (1/hq)*
          ((x[kn]-x[ks])/(2*hp))*
          (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
	        +
          (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
      insert(R, 1, &v, &C);
                
      C = ke;
      v = (0.5/(hq*hq))*
          (x[ke]-x[kw])*
          ((x[kn]-2*x[k]+x[ks])/(hp*hp))
          -
          (1/hq)*
          ((x[kn]-x[ks])/(2*hp))*
          (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
          +
          (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
      insert(R, 1, &v, &C);
                
      C = kn;
      v = (1/(hp*hp))*
          (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
          -
          (1/hp)*
          (x[ke]-x[kw])/(2*hq)*
          (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
          +
          (0.5/(hp*hp))*
          (x[kn]-x[ks])*
          (x[ke]-2*x[k]+x[kw])/(hq*hq)
          +
          Gamma(hp*r)*
          (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
      insert(R, 1, &v, &C);


      C = ks; 
      v = (1/(hp*hp))*
          (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
          +
          (1/hp)*
          (x[ke]-x[kw])/(2*hq)*
          (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
          -
          (0.5/(hp*hp))*
          (x[kn]-x[ks])*
          (x[ke]-2*x[k]+x[kw])/(hq*hq)
          -
          Gamma(hp*r)*
          (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
      insert(R, 1, &v, &C);
                                
      C = knw;
      v = (0.5/(hp*hq))*
          (x[kn]-x[ks])/(2*hp)*
          (x[ke]-x[kw])/(2*hq);
      insert(R, 1, &v, &C);

      C = kse;
      insert(R, 1, &v, &C);

      C = kne;
      v = -v;
      insert(R, 1, &v, &C);
      
      C = ksw;
      insert(R, 1, &v, &C);
    }
  }


  // derivative of rightmost column
  for (r=1; r<M-1; ++r)
  {
    k   = (r+1)*N  -  1;
    kw  = k-1;
    kn  = r*N      -  1;
    knw = kn - 1;
    ks  = (r+2)*N  -  1;
    ksw = ks - 1;
    ke  = r*N;
    kne = (r-1)*N;
    kse = (r+1)*N;
  

    R = C = k;
    v = (-2/(hp*hp))*
        (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
        -
        (2/(hq*hq)*
        pow((x[kn]-x[ks])/(2*hp), 2.0));
    insert(R, 1, &v, &C);
    
    C = kw;
    v = (-0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k]+x[ks])/(hp*hp))
      +
      (1/hq)*
      ((x[kn]-x[ks])/(2*hp))*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = ke;
    v = (0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k]+x[ks])/(hp*hp))
      -
      (1/hq)*
      ((x[kn]-x[ks])/(2*hp))*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn]-x[ks])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = kn;
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      +
      (0.5/(hp*hp))*
      (x[kn]-x[ks])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      +
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
    insert(R, 1, &v, &C);
    
    C = ks; 
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      +
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw]+x[ksw]-x[kse])/(4*hq*hp)
      -
      (0.5/(hp*hp))*
      (x[kn]-x[ks])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      -
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn]-x[ks], 2.0);
    insert(R, 1, &v, &C);
    
    C = knw;
    v = (0.5/(hp*hq))*
      (x[kn]-x[ks])/(2*hp)*
      (x[ke]-x[kw])/(2*hq);
    insert(R, 1, &v, &C);
    
    C = kse;
    insert(R, 1, &v, &C);
    
    C = kne;
    v = -v;
    insert(R, 1, &v, &C);
      
    C = ksw;
    insert(R, 1, &v, &C);
    
    insert(R, 1, &v, &C);
  }



  // derivative of bottom row
  for (c=1; c<N-1; ++c)
  {
    k   = (M-1)*N   +  c;
    kn  = (M-2)*N   +  c;
    ke  = k+1;
    kw  = k-1;
    knw = kn+1;
    kne = kn-1;


    R = C = k;
    v = (-2/(hp*hp))*
        (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
        -
        (2/(hq*hq)*
        pow((x[kn])/(2*hp), 2.0));
    insert(R, 1, &v, &C);
    
    C = kw;
    v = (-0.5/(hq*hq))*
        (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      +
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = ke;
    v = (0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      -
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);
    
    C = kn;
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (0.5/(hp*hp))*
      (x[kn])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      +
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn], 2.0);
    insert(R, 1, &v, &C);
    
    C = knw;
    v = (0.5/(hp*hq))*
      (x[kn])/(2*hp)*
      (x[ke]-x[kw])/(2*hq);
    insert(R, 1, &v, &C);
     
    C = kne;
    v = -v;
    insert(R, 1, &v, &C);
  }



  // derivative of bottom-left corner 
  {
    k   = (M-1)*N;
    kn  = (M-2)*N;
    ke  = k + 1;
    kw  = k + N-1;
    kne = kn + 1;
    knw = k - 1;

    R = C = k;
    v = (-2/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (2/(hq*hq)*
       pow((x[kn])/(2*hp), 2.0));
    insert(R, 1, &v, &C);

    C = kw;
    v = (-0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      +
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);

    C = ke;
    v = (0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      -
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);

    C = kn;
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (0.5/(hp*hp))*
      (x[kn])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      +
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn], 2.0);
    insert(R, 1, &v, &C);

    C = knw;
    v = (0.5/(hp*hq))*
      (x[kn])/(2*hp)*
      (x[ke]-x[kw])/(2*hq);
    insert(R, 1, &v, &C);

    C = kne;
    v = -v;
    insert(R, 1, &v, &C);
  }



  // derivative of bottom-right corner
  {
    k   = M*N - 1;
    kw  = k-1;
    ke  = (M-1)*N;
    kn  = ke - 1;
    knw = kn - 1;
    kne = (M-2)*N;

    R = C = k;
    v = (-2/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (2/(hq*hq)*
       pow((x[kn])/(2*hp), 2.0));
    insert(R, 1, &v, &C);

    C = kw;
    v = (-0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      +
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);

    C = ke;
    v = (0.5/(hq*hq))*
      (x[ke]-x[kw])*
      ((x[kn]-2*x[k])/(hp*hp))
      -
      (1/hq)*
      ((x[kn])/(2*hp))*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (1/(hq*hq))*(pow((x[kn])/(2*hp),2.0));
    insert(R, 1, &v, &C);

    C = kn;
    v = (1/(hp*hp))*
      (1 + pow((x[ke]-x[kw])/(2*hq), 2.0))
      -
      (1/hp)*
      (x[ke]-x[kw])/(2*hq)*
      (x[kne]-x[knw])/(4*hq*hp)
      +
      (0.5/(hp*hp))*
      (x[kn])*
      (x[ke]-2*x[k]+x[kw])/(hq*hq)
      +
      Gamma(hp*r)*
      (3.0/(8.0*pow(hp,3)))*pow(x[kn], 2.0);
    insert(R, 1, &v, &C);

    C = knw;
    v = (0.5/(hp*hq))*
      (x[kn])/(2*hp)*
      (x[ke]-x[kw])/(2*hq);
    insert(R, 1, &v, &C);

    C = kne;
    v = -v;
    insert(R, 1, &v, &C);
  }

  firstFill = false;
  J->FillComplete();
  return true;
}
